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Abstract 



We investigate the phase-space for trajectories in multi-black hole spacetimes. 
We find that complete, chaotic geodesies are well described by Lyapunov 
exponents, and that the attractor basin boundary scales as a fractal in a 
diffeomorphism invariant manner. 
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Chaos is an aspect of General Relativity which has only recently been explored. Its 
existence is expected as General Relativity contains generalizations of well-known chaotic 
systems from Newtonian physics such as the 3-body problem. Moreover, the non-linearity 
of Einstein's equations may give rise to chaos in systems whose Newtonian analogue is non- 
chaotic. The characterization of chaos in General Relativity is complicated by the dynamical 
nature of time, which necessitates a careful generalization of the parameters which quantify 
chaos, such as Lyapunov exponents. 

Relativistic systems in which chaos is known include charged particles in a magnetic 
field interacting with gravitational waves JTJ and particles near a magnetized black hole ||, 
as well as particles in Majumdar-Papapetrou (MP) geometries, described below. There is 
also the general, sufficient but not necessary, criterion that chaotic geodesies will occur if 
the phase space is compact and the Ricci scalar is negative || . In MP spacetimes the Ricci 
scalar is identically zero, so this criterion is inapplicable. 

Chaos is also found in cosmological models such as the Mixmaster universe 0, 0, 
and Robert son- Walker models containing dynamical fields JFj, however there is difficulty in 
defining coordinate invariant measures of chaos for such systems ||. 

MP geometries 0, |1(J are the static solution of the Einstein-Maxwell equations with 



metric and electrostatic potential given by 

ds 2 = -U- 2 dt 2 + U 2 (dx 2 + dy 2 + dz 2 ) , (1) 

A t = U' 1 , (2) 

and the spatial components of are zero, with U a function of the spatial coordinates 
satisfying Laplace's equation 

V 2 U (x, y, z) = U xx + U yy + U zz = . (3) 



Hartle and Hawking [[Tl|] showed that if U is of the form 

N M- 

t/ = l + E / (4) 
i= i \J{x- Xi) 2 + {y- y.j) 2 + (z- Zi) 2 

the MP metric corresponds to a system of extreme Reissner-Nordstrom black holes with 
equal charge and mass Mj > and horizons at (xi,yi,Zi). Note that the MP coordinates 
are singular at these points, mapping a horizon of finite proper area to a single point. They 
also showed that if U takes a form different to the above expression, the space-time contains 
naked singularities. 

Chandrasekhar |CJ and Contopoulos [I3|, have investigated the timelike and null 



geodesies of the two black hole system from the point of view of the periodic orbits and 
the weak field limit. For particles with elliptic (bound) energies the trajectories fall into 
several categories. There are stable periodic and quasiperiodic orbits, chaotic orbits trapped 
between periodic orbits, trajectories which fall into one or other of the black holes, and 
chaotic trajectories which lie on the boundary of these regions. This structure, together 
with the fact that the weak field limit is integrable, makes the MP geodesies problem a 
particularly interesting example of chaos in General Relativity. 
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Working from the equations of motion, Contopoulos showed that the weak field limit 
of the two black hole system was integrable for uncharged test particles with zero angu- 
lar momentum about the axis of symmetry. Using Hamilton- Jacobi methods ||15||, we have 



generalized this result to include test particles of arbitrary charge, energy and angular mo- 
mentum. In the background spacetime generated by two black holes with masses Mi, M 2 
centered at (0, 0, 0) and (0, 7r, 0) in prolate spheroidal coordinates, (ip, 9, <fi) [|13|], we find the 
Super-Hamiltonian for a test particle of charge e and mass m is given by 



H 



dS_ 



'dS_ _ _e x 
~dt~U 



2U 2 Q 




1 fdS_\' 
2(f/sinh^sin^) 2 J 



(5) 



where 



U = l + W/Q , Q = sinh 2 ^ + sin 2 Q , 
W = (Mi + M 2 ) cosh rj> + (Mi - M 2 ) cos 6 . 

Since t and are cyclic coordinates, the canonically conjugate momenta n t = —{E — e), 
and 7r^ = L z are constants of the motion. E and L z are the particle's energy and angular 
momentum at infinity. Invariance of the particle's rest mass along the trajectory aflinely 
parameterized by A gives rise to a third constant of motion, m 2 . We see that the system is 
separable, and hence integrable, if 



S 



(6) 



(7) 



Substituting this expression into (|) and using the weak field expansion 

QU n = Q + nW + 0(W 2 /Q) , 

we find 

= J (^(E 2 + m 2 ) sinh 2 ip + 2(M 1 + M 2 )(2E 2 + m 2 - Ee) cosh^ - 

S = J (^(E 2 + m 2 ) sin 2 9 + 2 (Ml - M 2 )(2E 2 + m 2 - Ee) cos 6 - 

The separation constant a is a fourth constant of the motion, which in turn guarantees the 
system is integrable. 

The next order in the weak field expansion destroys separability in prolate spheroidal 
coordinates. While this does not prove the system is non-integrable, it does show that if 
chaos is present, it is due to relativistic effects. In this case, the chaos is related to the 
relativistic perihelion advance of elliptic trajectories, which invalidates Bonnet's theorem for 




motion with two centres 12 



For MP geometries with three or more black holes, even the Newtonian limit is non- 
integrable. The one exception to this occurs in the Newtonian regime for "extremal" test 
particles with e = m, since then the gravitational and electrostatic forces cancel exactly. 
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This cancellation breaks down in General Relativity where gravity couples to all forms of 
energy, including kinetic energy. The scattering of an extremal test particle by an extremal 
black hole of mass M illustrates this effect since we find the scattering angle A0 differs 
from the Newtonian value (A0 = 0) by A0 = -(M/6)[3 + u 2 /A + 0{u A , M/b)]. The impact 
parameter b is related to L z and the asymptotic velocity u via L z = mbuj \J\ — u 2 . This 
result leads us to expect that the trajectories of extremal test particles in MP geometries 
will be chaotic beyond the Newtonian limit, thus providing another example where chaos 
arises due to relativistic effects. 

For rapid numerical calculations we use the original MP coordinates (t,x,y,z), because 
the equations of motion do not contain transcendental functions. Writing the upper com- 
ponents of the four- velocity in an orthonormal basis as (t>o,v), the equations of motion 
are 

ir = U~ 2 [(vl + v 2 -ev /m)VU -vv ■ W] , (10) 
i=U~ 1 v , (11) 
t = U vp , (12) 
v = Vl + v 2 , (13) 

where a dot indicates the derivative with respect to proper time. The conserved energy is 
given by 

E = U- 1 (mv -e) . (14) 




FIG. 1. A two-dimensional section of phase space determined by v = 0, y = 0. The black 
regions are trajectories which fall into the black hole of mass 1/3 at (0, 1). The white regions are 
those which fall into the black hole of mass 1/3 at (0, —1). 

We have integrated the equations using a 4th order Runge-Kutta technique with adaptive 
stepsize |16|]. The results are shown in Figs. 1-3. These figures show two-dimensional sections 
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of the phase space given by the initial conditions v = 0, y = 0. All of the points shown 
fall into one of the black holes and are color coded accordingly. The boundary between the 
basins of attraction clearly has a complicated structure, which may be quantified using the 
concept of fractal dimension. 
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FIG. 2. A small region of Fig. 1 showing the fractal structure. 

There are a number of different fractal dimensions defined in the literature ||17|| . The 
box-counting dimension of a set embedded in an .E-dimensional Euclidean space is obtained 
by finding the minimum number N(e) of i?-dimensional cubes needed to cover the set, and 
calculating 

In N 

d B = -lim-j , (15) 

e-»o ine 

assuming that this limit exists. The generalization to Riemannian or pseudo-Riemannian 
spaces is direct, since d B is invariant under diffeomorphisms, so all choices of coordinates 
give the same value. If the basin is defined on some spatial hypersurface, as is the case here, 
the dimension may depend on the chosen slicing of spacetime. 

The value of d B of the fractal boundary for the region of phase space shown in Fig. 2 
was evaluated using the above equation. The region, containing 2520 2 points, was covered 
by a grid. Each square, of size e 2 (e being a factor of 2520) was counted if it contained 
trajectories leading to both black holes. Fig. 4 shows a plot of IniV vs lne. The straight 
line is a fit to all but the three largest and smallest values of e, and gives a dimension of 
1.464 . . .. This result clearly demonstrates that the 8-dimensional boundary as a whole is a 
fractal, thus giving us a diffeomorphism invariant measure of chaos. 
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FIG. 3. The v = 0, y = section of phase space for the 3-hole problem. The black holes of 
mass 1/3 are situated at (0,1), (y/3/2, — 1/2), and (— V3/2, — 1/2), corresponding to gray, black, 
and white, respectively. 
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FIG. 4. The dimension ds of a region of the fractal boundary near that shown in Fig. 2 is 
equal to 1.464 .... 



For complete geodesies, a useful characteristic of chaotic systems is the presence of pos- 

These 



itive Lyapunov exponents, which demonstrate that the system is nonintegrable [18 



are defined in flat space-time by choosing a point x in phase space, at the centre of a ball 
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of radius e C 1. After a time t the ball evolves to an ellipsoid with semi-axes where k 
ranges from one to the dimension of the phase space. The Lyapunov exponents are 



\k(x) = lim lim - In k , (16) 

again assuming the limits exist. The are constant along a trajectory, and are often 
constant over larger regions of phase space, such as the basin of an attractor. Numerically 
these are calculated by integrating the linearized equations and periodically performing a 
Gram-Schmidt orthonormalization to ensure that the exponentially growing solutions do not 



cause overflow errors [19] 



There are a number of subtleties associated with the definition and calculation of Lya- 
punov exponents in curved spaces, and in particular, in the MP space-times. The definition 
of the Afc requires a metric on phase space to calculate distances, but the Xk are independent 
of this metric if the phase space is compact, since the trajectory passes arbitrarily close to 
at least one point an infinite number of times at unbounded values of t, causing the metric 
terms appearing in the ratio of radii to effectively cancel out. 

In relativity, there is no unambiguously defined global time in general. If this is the 
case, the only reasonable option is to use the particle's proper time as t in Eq. flnS|). This 
approach emphasizes the quasi-local nature of Lyapunov exponents, that is, their relation 
to particular trajectories, and not (necessarily) the system as a whole. If the system has a 
timelike Killing vector, and hence a global time t, then this may be used in the definition 
of the exponents. The results of this approach correspond to measurements of stationary 
observers. We choose the latter approach, following Ref. [0]. 

In the MP case, there is an additional difficulty in that many of the geodesies are incom- 
plete. Particles may cross the horizon of a black hole after finite proper time and encounter 



a singularity, rendering the limit in Eq. (|16"D undefined. Thus, we attempt to calculate 



Lyapunov exponents only for those trajectories which survive for an infinite proper time. 
This is not a problem for the regions of phase space which are nearly integrable, but for 
the basin boundary special techniques are required. It is impossible to find a point which is 
exactly on the boundary, and a trajectory which begins near the boundary will eventually 
move away and fall into a black hole. The solution is to periodically adjust the trajectory 
to ensure that it remains close to the boundary. This is achieved by small random shifts, 
checking to see that the resulting trajectory survives more than a specified proper time. We 
have checked that this procedure is stable with respect to the precision used and does not 
result in systematic shifts in the energy E over the integration times used. We have calcu- 
lated the Lyapunov exponents for the two-black hole system using the 4-dimensional sub- 
manifold of phase space, (x, z, v x , v z ), defined by [L z = 0, m = const), for the initial values 
(3.33467,0.23509,0,0), and found them to be \ k = (0.03609,0.00006,-0.00006,-0.03609). 
The statistical uncertainty is about 0.00009. Note that the Lyapunov exponents are mea- 
sured in the eigendirections of the linearized evolution operator, and do not correspond to 
the above coordinate system. 

As in non-relativistic mechanics, the symmetries of the system impose constraints on the 
Lyapunov exponents [ 2"0"fl . Liouville's theorem implies that they sum to zero; time reversal 



symmetry implies that they come in +/— pairs; and the presence of a constant of the motion 
requires that one pair of the exponents be zero. Since the phase space for our system is eight 
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dimensional, we require four constants of motion to ensure the vanishing of all four pairs 
of Lyapunov exponents. For geometries with three or more black holes, there will generally 
only be two constants of motion, necessitating the calculation of two pairs of Lyapunov 
exponents. These results will be presented in a future study. The two-black hole geometry 
furnishes three constants of motion, leaving one pair of non-zero Lyapunov exponents as 
given above. 

Finally, we note that recently the MP solutions have been generalized to include a positive 
cosmological constant pif , |2^| . These solutions may describe coalescing extremal black 



holes in a de Sitter type universe. It would be interesting to investigate this system, and 
how its dynamical nature and negative Ricci curvature might affect the chaotic structure of 
the geodesies. 
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